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Lattice models of coupled dynamical systems lead to a variety of complex behaviors. 
Between the individual motion of independent units and the collective behavior of mem- 
bers of a population evolving synchronously, there exist more complicated attractors. In 
some cases, these states are identified with self-organized critical phenomena. In other 
situations, with clusterization or phase-locking. The conditions leading to such different 
behaviors in models of integrate-and-fire oscillators and stick-slip processes are reviewed. 
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Criticality. 

1. Introduction 

In nature, there are many examples of systems with complex collective behavior. 
There is a big effort to identify and understand the underlying mechanisms leading 
to such behavior, but it is difficult to find general rules which could give, a priori, 
information about spatio-temporal complexity. The analysis of the time evolution of 
the magnitudes which describe an isolated dynamical system is a first step. However, 
the features of a big system consisting of a large number of individual units (where 
each unit is a dynamical system itself) interacting according to a given criterion, 
can be very difficult to determine, even if all the units are identical. The nature of 
the interaction between units, the type of boundary conditions, and the absence or 
existence of noise are some ingredients which can change completely the dynamic 
behavior of a coupled dynamical system. 
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One might think that complexity arises when the intrinsic dynamics that govern 
the temporal behavior of each member of a population as well as the interaction be- 
tween them follow complicated spatio-temporal rules. Not necessarily. It may arise 
as a result of the continued local simple interactions between all parts in an extended 
system. On the other hand, extended systems with complex local features may lead 
to a "simple" collective behavior. This is observed in some biological systems where 
after some transient period a regime characterized by a perfect synchrony in the 
temporal activity of all the members is achieved. 

Between the trivial behavior of independent uncorrelated units, which can be 
analyzed by considering the features of isolated elements, and "simple" collective 
behaviors such as units evolving in perfect synchrony, there is a wide spectrum of 
situations. Some of them manifest a large degree of complexity such as systems dis- 
playing spatial fractal structures and temporal fractal behavior. Fractal structures 
appear in a variety of physical systems. Era They are geometrical objects, that look 
alike on all length scales. The universe consists of clusters of galaxies, organized in 
clusters of clusters of galaxies and so on. Earthquakes occur on structures of faults 
ranging from thousands of kilometers to millimeters. Likewise, the coast of Norway 
contains fjords of all sizes, from very small ones to very big ones. Temporal fractal 
behavior occurs when the temporal fluctuations in a certain quantity look alike on 
all time scales. In general, the power spectrum behaves like 1//^. When ip w 1 
we talk about 1/f noise or flicker noise.QoO The fluctuations in light intensity of 
quasars, the flow of the river Nile, and the current flowing through a resistor are but 
a few examples of systems displaying fluctuations without an intrinsic time scale. 

During the last years there has been an increasing interest to understand the 
mechanisms that lead to other complex behaviors closely related to fractals. In 
particular, a lot of attention is paid on conservative and nonconservative systems 
displaying self-organized criticality (SOC). Up to now, SOC has been observed in 
models of cellular automata, coupled map lattices, and coupled dynamical systems 
defined on a lattice. It is our purpose in this paper to review the mechanisms as 
well as the conditions required to observe SOC and the corresponding transition to 
other collective behaviors. 

The paper is organized as follows. In Section II we review the mechanisms that 
provoke the occurrence of SOC, in randomly driven as well as in deterministically 
driven models. In Section III we analyze different models of oscillators focusing our 
interest on the conditions required to observe a macroscopic degree of synchroniza- 
tion. We will discuss the effect of both short- and long-range interactions. Section 
IV is devoted to investigate the circumstances in which a system develops the be- 
haviors mentioned above, their transitions, and their relations. Finally, in Sect. V 
we present the conclusions of the present work along with the lines we think deserve 
further investigation in order to deal with more realistic physical problems. 

2. Self-Organized Criticality 

The ubiquity of spatial fractal structures and temporal fractal behavior observed 



Self- Organized Criticality and Synchronization in Lattice Models 3 



in physical systems suggests a common underlying dynamical origin. The lack of 
a characteristic scale is the hallmark of a critical process. In equilibrium, critical 
processes require a fine-tuning of a relevant physical control parameter such as the 
temperature or the magnetic field to a critical point. However, nature does not, by 
itself, provide any fine-tuning of control parameters, and with zero probability it sits 
at the critical point by accident. Thus it is very unlikely that the wide occurrence 
of scale-invariance is due to critical processes in equilibrium systems. 

Recently, Bak, Tang, and Wicscnfcld (BTW) suggested that the frequently ob- 
served scale-invariance in nature might be related to the spontaneous organization 
(spatial complexity and reorganization (temporal complexity) in slowly driven, dis- 
sipative systems.Bi3 Internal mechanisms lead to dissipation (heat generation) and 
quasi-stationary states are reached in which the average fluxes of energy into and 
out of the systems are equal. The avalanches that occur when grains are dropped 
onto a pile have been used to illustrate this idea. When sprinkling grains of sand 
onto a table, one after the other, a pile builds up. Eventually, the pile ceases to 
grow and additional grains of sand will ultimately fall off the pile. The attractor of 
the dynamics is a statistically stationary state with a fluctuating angle of repose. 
If the slope is too steep, large system spanning avalanches would make the pile 
to collapse to a more stable configuration. If, on the other hand, the slope is too 
shallow, only small avalanches would be initiated. The sandpile settles into a state 
in between with avalanches of all sizes, power-law distributed. This scale invariancc 
suggests that the system is critical in analogy with equilibrium critical phenomena. 
However, one deals with dynamical nonequilibrium statistical properties and the 
system evolves naturally to the critical state without any tuning of external control 
parameters. The criticality is an intrinsic property of the dynamics of the system, 
and the phenomenon of self-organized criticality may very well provide a connec- 
tion between the occurrence of fractal structures and 1// noise, as well as being 
the physical origin of these two phenomena. The earth crust is another example 
of a self-organized critical system: Along the boundary of tectonic plates, a slow 
build-up of strain is relaxed through earthquakes of all sizes; the energy-frequency 
relation of earthquake occurrence, which is related to the Gutenberg-Richter lawja 
is a power-law distribution. 

A cellu l a r automaton model for sandpile dynamics&B and other numerical 
modelsoii30H3 were shown numerically to display SOC. The dynamical rules in 
the BTW model — also known as the sandpile model — at least intuitively re- 
semble the dynamics of a sandpile. In the sandpile model, each site on a lattice 
is characterized by an integer variable (local slope), and this variable changes in 
time due to external perturbations (adding grains of sand) and to interaction be- 
tween different sites when an avalanche propagates through the system. With very 
simple rules, the system reaches a state which is characterized by the lack of any 
characteristic length or time scale. The simplicity of the model suggests that the 
phenomenon of SOC could be quite universal: Indeed, it has_been scrutinized by 
a host of researchers in fields spanning statistical mechanics, tij condensed matter 
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The notion of self-organized criticality has initiated much experimental work, 
especially on granular systems. The main object has been to address the question 
of SOC in real sandpiles. Essentially, two distinct types of experiments have been 
performed: Rota tin g a s e n - ii-cylindrical drum partially filled with grains at a low, 
constant velocityj!3l2j'El[E3 or droppine^singlegrains at a low rate on a conical 
sandpile resting on a circular platform.r 3 H 24 iT 5 tr 6 l In most of these studies, the flow 
over the rim of the system (the drop number) was recorded. In a drum experiment, 
this quantity was found to be nearly periodic in time.113 Similar behavior was seen 
for large conical piles: they displayed relaxation oscillations when a system size of 
30-40 particle diameters was exceeded.LJ In some of the latter experiments, power- 
law distributions has been reported for small conical piles. However, recently it has 
been shown, that these distributions are more likely to be stretched exponentials, 
that is, a characteristic size for the drop number appears. E3c3 The observed behavior 
may be attribute d t o iner tia! effects, which seem to play an important role in all 
the experiments.E^OcSEj Inertial effects lead to a nonlocal process, whereas in the 
numerical models only the local geometry determines the dynamics. Furthermore, 
we notice, that the relationship between the drop number measured in experiments 
and the avalanche size measured in numerical models is unclear. For an alternative 
view on the dynamics of granular material see the work by Mehta and Barkertil 
which reviews the progress that has been achieved experimentally and theoretically. 

In a recent experiment, the internal dissipated energy_due to avalanches in a 
slowly driven one-dimensional (ID) rice pile, was recordedcH and it was shown that 
the occurrence of SOC depends on details in the grain-level dissipation mechanism. 
With spherical grains, a stretched-exponential distribution was observed implying 
a characteristic scale which is inconsistent with SOC. The spherical grains typically 
accumulated kinetic energy when moving down the slope. However, with more 
elongated grains, the dynamics was dominated by sliding grains. This induced a 
higher effective friction which suppressed the inertial effects, and a power-law dis- 
tribution of avalanche sizes appeared. This provides the first experimental evidence 
of self-organized critical behavior in slowly driven granular systems. 

2.1. Randomly driven models 

First, we briefly discuss the ID BTW sandpile model. Although the behavior of the 
ID model is trivial, it illustrates the basic concepts of the more complex behavior- 
displayed by the 2D BTW sandpile model. 

The ID sandpile modefl is a cellular automaton where an integer variable hi 
gives the height of the pile at site i. We define the local slope Zi at site i as 

Zi = hi-h i+ i. (1) 

The addition of a grain of sand on a randomly chosen site i (hi — > hi + 1) results in 
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the following changes in the slopes 

z ^ 7 ZV 1 ( 2 ) 

Zi-l —> z i-l — J-- 

We proceed by dropping grains at random sites until one site reaches a slope which 
exceeds a critical value, Zi > z c , and the site topples by transferring one grain to 
its neighboring site on the right, that is, 

Zi * %i 2 , , 

Zi±l ->■ Zi±i + 1 

unless at the rightmost site where the sand grains fall off the pile. The neighbors 
that are affected by the toppling may topple in turn generating a chain reaction or 
avalanche. During the avalanche, no grains are added to the pile thus separating 
the two time scales involved in the dynamic evolution of the pile, a slow time scale 
for the addition of the grains and a fast time scale for the relaxation processes. 
The avalanche stops when the system reaches a stable state with zi < z c Vj and 
another grain is added following (^) until a new avalanche is initiated and so on. 
After a transient period, whose duration depends on the initial conditions, the 
system reaches a critical state in which — z c Vi. This state is a fixed point 
of the dynamics, since after any perturbation the system relaxes returning to the 
stable state; the added grain will tumble down the slope and simply fall off the pile. 
The fixed point is an attractor for the dynamics, however, this state has no spatial 
structure, and correlation functions are trivial. 

In the 2D sandpile model (on a square lattice), El an integer Zjj is assigned 
to each lattice site where i,j — 1,...,L. The integer Zij represents an 

appropriate dynamical variable (e.g. height of a column of sand, mechanical stress, 
heat, pressure, or energy^El) on site in a spatially extended system. 

We perturb the system (add sand to it) by choosing at random a position and 
increasing the dynamical variable with one unit, that is, Zi : j — > Zij + 1. Whenever 
the dynamical variable on site exceeds a threshold value Zij > z c , the site 

topples, 

Zij ► Zij — 4, (4) 

Znn * Z nn -\- 1 (5) 

where z nn denotes the height at nearest neighboring sites. As a result one or more 
neighbors may exceed the threshold value in which case they have to relax and 
an avalanche will propagate through the system. The toppling rule conserves the 
amount of z-values whenever an interior site topples. Dissipation only occurs when 
a boundary site topples assuming open boundary conditions (be). 

The 2D sandpile will evolve into a statistically stationary state where, on the 
average, the rate of flow into the system equals the rate of flow out of the system 
across the boundary. Whenever the average slope becomes too large, system span- 
ning avalanches can occur, which transfer grains to the boundary. On the other 
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hand, when the average slope is too small, avalanches tend to be small, and the pile 
builds up. The stationary state is no longer a fixed point as in ID but a complex 
attractor. DharE3 has calculated the exact number of states in the attractor. When 
z c = 3, only (3.210 • • -) L of the 4 L stable states are allowed. Furthermore, the 
dynamics which takes the system from one allowed state to another is ergodic in 
the sense that all the allowed states occur with equal probabilityEj 

One way of characterizing the dynamics in a system of linear size L, is to measure 
the distribution of avalanche sizes defined as the total number of topplings. Let 
P{s, L) denote the probability of initiating an avalanche of size s in a system of size 
L. In the stationary state, an added particle will eventually fall off the boundary 
of the system. The average distance a particle has to flow in order to reach the 
boundary is proportional to L due to the random deposit. The average avalanche 
size (s) = J sP(s, L) ds must be infinite in the thermodynamic limit L — > oo. 
Indeed, averaged over a large number of perturbations (after the transient period), 
the avalanche size distribution is a power law 

P(s 1 L)<xs 1 - T (6) 

limited only by the size of the system, see Fig. |l](a) . From these results the power- 
law exponent r is 2.15 ± 0.1. Since the cutoff is a finite size effect, the average 
avalanche size tends to infinity when L — > oo. A data collapse for different system 
sizes L is obtained when plotting LP P(s, L) or s T ~ 1 P(s, L) against the rescaled 
variable s/L v , see Figs. |l|(b-c). Thus we can write 

P{s,L)=L-0F{s/L v ), (7) 

or, alternatively, 

P{ S ,L) = s 1 - T G{ S /L») 1 (8) 

with G(x) = x T ~ 1 F(x). The exponent (3 = v(t — 1), where r is the power-law 
exponent and v a critical index expressing how the cutoff scales with system size. 
The scaling function F approaches a power law (or, equivalently, G a constant) 
when s/L v — > since the avalanche size distribution becomes independent of the 
system size when L — > oo and decays very quickly for s» L u . 

The BTW relaxation rules (^) and (^) conserve the dynamical variable z except 
at the boundary. Introducing nonconservative dynamics in the interior of the system 
will leave the BTW model noncritical. Changing rule (|]) to, say, Zjj — » Zi t j — 5 
while (§) remains unchanged, one grain of sand will dissipate for every toppling. In 
the stationary state, the average rate of dissipation equals the average rate of flow 
into the system. Thus (s) = l/P Za where P Za is the probability of adding a grain to 
a site with z c units, thereby initiating an avalanche. This probability will approach 
a constant value when L — > oo, that is, (s) will approach a finite value. Such a 
system cannot display SOC. The avalanche size distribution decays exponentially 
with a characteristic avalanche size, that does not depend on the system size, and 
the avalanches are "localized". However, a system with globally conservative but 
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Fig. 1. The 2D Bak, Tang, and Wiesenfeld sandpile model for system sizes L = 50, 100, and 
200 with open be. (a) The probability of initiating an avalanche of size s in a system of size L, 
P(s, L) decreases algebraically with s. The cutoff is a finite-size effect. The displayed distribution 
functions are averaged over exponentially increasing bins, (b) Using Eq. (ft, a reasonable data 
collapse is obtained with r = 2.15 ± 0.1 and v = 2.7 ± 0.1. (c) However, the alternative data 
collapse suggested by Eq. (H) shows the bad scaling of the model for "small" system sizes. 
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Fig. 1. (continued) 

locally nonconservative dynamical rules (changing rule (^) to Zij — ► Zij —4 + 9, 
where 9 G {—3, 3} is an annealed random variable) displays SOC.0 

Memory effects can be introduced into the BTW model in the sense that a site 
which has relaxed cannot receive any amount from neighboring sites during the on- 
going avalanche. This kind of memory effects makes the dynamics nonconservative 
and hence destroys the SOC behavior. 

In order to model the early experiments on SOC in slowly driven granu- 
lar piles, inertia effects were included in a BTW-like modelEJ It was observed 
that when the system size exceeded a certain size, characteristic system spanning 
avalanches ap peared du e to the inertia effects in agreement with experiments on 
real sandpiles.llaEj'Ea'EJ Thus the presence of inertia effects also destroy the SOC 
behavior. 

An apparently different stochastically driven continuous-energy model was in- 
troduced by Zhang.E3 In this model, a real number £!y is assigned to each site (i, j) 
on a square lattice. The system is perturbed randomly by adding an amount SE to 
a randomly chosen site, that is, Eij — > Eij + SE. If the value at any site exceeds 
the threshold value Eij > E c , the system relaxes according to 

Enn * E nn -\- cEi^j , * 

with e = 0.25. The system is conservative and dissipation occurs only at the 
boundary, where the number of nearest neighbors is less than 4. The model was 
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studied with E c = 1 and SE chosen uniformly in the interval [0, 0.5],l3E3 in which 
case it is believed to be in the same universality class as the 2D BTW model, that 
is, the values of critical exponents, such as the power-law exponent t and scaling 
index v, are identical.EjO 

With e < 0.25, the dynamics in the Zhang model is nonconservative. When 
SE is chosen in the range [0.1, 1], the nonconservative relaxation rule introduces a 
characteristic avalanche size. In this case, the Zhang model is also noncritical in the 
presence of nonconservation. Furthermore, the introduction of the memory effects 
mentioned above leads to a characteristic avalanche size inconsistent with the hy- 
pothesis of SOC. Similar models have been considered in the context of fracturingHa 
However, additional rules where introduced in order to make the dynamics globally 
conservative which introduces a kind of inertia effects. Indeed, the general behavior 
of such systems is the periodic generation of system spanning avalanches coexisting 
with power-law distributed small avalanches. 

The microscopic rules of the lattice models can be written in the form of a set of 
coupled equations, one for each site. For a continuous version of the stochastically 
driven BTW model we haveEa 

E id (t + 1) = E id (t) - e(E itj (t) - E c ) E c + Y,®{E nn {t) - E c ) E c /A + rnj(t), (10) 

nn 

where 0(x) = when x < and O(x) = 1 when x > is the Heaviside function. 
The variable t refers to lattice updates. However, in a slowly driven system, the 
physical time between lattice updates is vera_short when avalanches are propagating 
and very long when perturbing the systemic The sum takes into account a possible 
increase in the dynamical variable due to toppling of nearest neighboring sites nn. 
When an avalanche has come to a halt, a random site is perturbed. The random 
drive is represented by 

mj(t) =Ti5 itn 5 j , m '[[[i - e(E h i(t) - e c )} , (ii) 

k,l 

n, m being two random integers between 1 and L chosen once for each lattice up- 
date, and k, I two index running over all the lattice sites. In the original BTW-model 
n = 1. The Heaviside function that appears in the noise terms represents the sepa- 
ration of time scales. The models are driven slowly in the sense that no perturbation 
takes place while avalanches evolve. On the other hand, for the conservative Zhang 
model, where the toppling variable is reset to zero, the corresponding equations are 

E id (t+1) = [l-e(E itj (t)-E c )]E itj {t)+J2Q(Enn{t)-E c )E nn (t)/4+r] itj (t). (12) 

nn 

The noise T}i,j(t) is again given by (pT|). In the original Zhang model, n is chosen 
uniformly in the interval [0,0.5]. In order to deal analytically with these equations, 
one usually neglects the Heaviside function in the noise term. This fact creates 
some difficulties, since the notion of avalanches in such models is not well-defincdo 
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Equations @ and @) can be coarse-grained in order to obtain a continuum 
equation for the effective E(r, t). Using the prescriptions for the temporal derivative 
and the Laplace operator one gets, after a rescaling of the energy E — E c 
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aV z [{ZE{r, t) + E c )e(E{r, t))} + r,(r, t), (13) 



where Z = 1 for the Zhang model and for the continuous BTW model. The noise 
7y(r, t) takes into account the effective external noise as well as the internal noise 
that appears due to the elimination of internal degrees of freedom. 

On the other hand the O-function representing the threshold dynamics is mod- 
eled as the limit of a smooth function, whereupon dynamical renormalization group 
calculations can be applied E3o This fact also changes the problem intrinsically. 
Left to itself, without adding noise, the model would relax to a unique flat ground 
state with E(r, t) = 0. The concept of stable states (or memory) is removed. Some 
other authors have built nonlinear stochastic differential equations according to the 
symmetries of the discrete models.Elo' H One should distinguish, however, between 
these approaches, in which the main goal is to study the "generic scale invariance" 
by means of dynamical renormalization group theory, and the lattice models. 



2.2. Continuously driven models 

A continuously driven dynamical system was introducecLhy Olami, Feder, and Chris- 
tensen (OFC) in the context of earthquake modelingHl Though the dynamics of 
earthquakes is very complex there are two basic components which have to enter a 
model: (a) earthquakes are generated by the very slow relative motion of tectonic 
plates, (b) they occur as abrupt rupture events when the fault can no longer sus- 
tain the stress, that is, the occurrence is intermittent. Hence, there are two time 
scales involved in the process; one is related to the stress accumulation while the 
other, which is orders of magnitude smaller, is associated to the duration of the 
abrupt releases of stress. A simplified spring-block model introduced by Burridge 
and KnopofEJ includes the basic components mentioned above. The model consists 
of a 2D network of blocks interconnected by springs. (For a review on the work of 
ID Burridge-Knopoff models of earthquake faults, see the work by Carlson, Langer, 
and Shaw £3). Each block is connected to the four nearest neighbors. Additionally, 
each block is connected to a single rigid driving (tectonic) plate by another set of 
springs as well as connected frictionally to a fixed rigid (tectonic) plate, see Fig. ||. 

Strain is accumulated uniformly across the system as the rigid plates move with 
a constant relative velocity. When the strain on one of the blocks exceeds the static 
friction force, the block slips. The released stress is transferred to the neighboring 
blocks, which in turn may slip, and an earthquake can evolve. This simple picture 
can be mapped to a coupled map lattice model.cH A uniformly increasing stress is 
assigned to each lattice site When a block exceeds the threshold stress E c 

it slips. Simple arguments lead to the relaxation rules identical to Eq. (^) with 
e = K/ (4K + Kl), K and Kl denoting the spring constants in the model, see Fig. 
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[|. The earthquakes are considered to be instantaneous, that is, the loading of the 
system stops during an earthquake and the driving recommences when the evolving 
earthquake has come to a halt. The size of an earthquake is defined as the total 
number of relaxations. 

When e = 0.25, the model is a continuously driven version of the Zhang model." 
Apparently, this model has the same power-law exponent r as the stochastically 
driven models. However, since Kl ^ 0, the generic situation corresponds to the 
nonconservative case where e < 0.25. 

In Fig. ||(a) we plot the results of a simulation with e = 0.20 for different system 
sizes. A reasonably good data collapse is obtained using the finite-size scaling ansatz 
Eq. (0), see Fig. ||(b). Thus the system is scale invariant; it displays SOC behavior 
even in the nonconservative case. 

When introducing nonconservative relaxation rules into the stochastically driven 
BTW and Zhang models, the distribution of avalanche sizes decay exponentially 
with a characteristic avalanche size. The systems are subcritical. Thus the occur- 
rence of criticality in the nonconservative OFC model is very intriguing since it 
suggests a different mechanism for the generation of the scale invariance. Since the 
majority of natural phenomena are inherently nonconservative, the SOC behavior 
of nonconservative systems is probably much more generic than the corresponding 
behavior of conservative systems. 



"There is also a uniformly driven continuous version of the 2D BTW sandpile modelfj Choosing 
random initial conditions, this model will be in the same universality class as the stochastically 
driven integer model. In the continuous model, the fractional part of the dynamical variable takes, 
in some sense, the role of the random number generator. Note that this model and the OFC model 
are deterministic. The inherent randomness enters via the initial conditions. 
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Fig. 3. (a) The probability density P(s,L) of having an earthquake of size s when e = 0.20, for 
the OFC model. The different curves refer to different system sizes L = 50, and 100. The cutoff 
in energy distribution increases with pyatq m size L. (b) A finite-size scaling plot using Eq. (|t]) 
with t = 2.9 ± 0.1 and v = 2.15 ± O.lffl 
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The microscopic rules for the continuously driven OFC model are given 




E itj (t + 1) = [1 - Q(E itj (t) - E c )\ Ei d (t) + eJ2 Q(E nn (t) - E c ) E nn (t) + v (t) (14) 

nn 

where 

v(t)=vl[[l-Q{E ktl (t)-E c )] (15) 

k,l 

is a global perturbation. In simulations, one uses rj = E c — E max (t) where E max (t) 
is the maximal value of all the dynamical variables in the lattice after the avalanche 
has come to a halt. To simulate this model a very efficient algorithm has been 
proposed by GrassbergerEJ 

The only difference between the Zhang model and the OFC model is the way of 
driving. The uniform (global) xlriving must be a crucial clement in driving noncon- 
servative models to criticalityH To move continuously between stochastically and 
uniformly driven systems, one can perturb the system randomly by adding a fixed 
amount SE to a randomly chosen site that is, E i: j — > E iy j+SE. If the value at 

any site exceeds the threshold value Eij > E c , the system relaxes according to Eq. 
([^). The limit of small SE is the limit of continuous drive. Figure ^ shows the dis- 
tribution function of avalanche sizes P(s, L) for e = 0.20 and various values of SE. 
The avalanches are localized for "large" values of SE. The distribution function is 
essentially an exponentially decreasing function and no scale invariance is observed 
when changing the system size. Apparently a phase transition from a localized into 
a critical system occurs between Fig. |(c) and |(d). For "small" values of SE the 
distribution function approaches a power-law behavior and the cutoff scales with 
system size, implying that the avalanches are no more localized. 

IE E (t) I 

The random drive represents an effective noise given by y — — ^ SE oc v SE. 
Obviously, large noise can destroy spatial correlations in a coupled dynamical sys- 
tem, b but t he S E — > represents the weak noise limit where long-range correlations 
can appear .til Analytically, the nonconservation introduces a length scale which in 
principle one expects to break the scale invariance. Additionally, the continuous 
driving is difficult to handle analytically. However, an appealing approach would 
be to treat the OFC model in the limit of random weak noise driving since this is 
conceptually equivalent to a continuous driving. 

3. Synchronization 

Synchronization is a fascinating phenomenon observed in biological, chemical, and 
physical systems that has captivated the attention of many scientists in the last 
century. The paradigmatic example is observed in some forests of South Asia where, 

6 This is seen for example in the model discussed by Manna, Kiss, and KerteszJO where the 
dynamical rule destroys the spatial correlations in the BTW model by letting the toppling site 
Zj j — » Zi j — 4 + 9, where 8 £ {—3, . . . , 3} is an annealed random variable. They obtain r s = 
2.515 ± 0.02|4fl|-«erfect agreement with the analytical result t s = 5/2 for models without spatial 
correlations .Mr 2 ! 




Fig. 4. The nonconscrvativc stochastically driven Zhang model with e = 0.20, E c = 1, and 
L = 25, 50 when (a) 5E = 1, (b) SE = 0.2, (c) SE = 0.008, and (d) SE = 0.00032. 
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Fig. 4. (continued) 
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at night, myriads of fireflies are sitting in the trees. At the beginning, the members 
of the population emit flashes of light incoherently but after a short period of 
time the whole swarm is flashing in unison creating one of the most striking visual 
effects ever seen. But mutual coherence in the temporal activity of a group of 
units has been observed in other contexts. Some examples are the cell mitotic 
cycle, heart pacemaker cells, circadian pacemakers (biological clocks), neurons in 
the visual cortex, menstrual period of women. Josephson junctions, or some chemical 
reactions, just to cite a few. Winfree's bookH is an excellent reference to read about 
these systems. 

The relevance of synchronization has been stressed frequently although not al- 
ways fully understood In the case of fireflies it may help the courtship between male 
and femalefir 5 ! In the heart, the impulses coming through the vagus nerve trigger 
the contraction of the heart only if they are properly synchronized £3 In other cases 
the relevance is still a matter of discussion. There are some evidences which support 
that coherent oscillations play an important role in sensory processing. It has been 
suggested that the discrimination and segmentation of objects, so often performed 
by living systems, can be explained hy a nalyz ing the temporal firing patterns of 



different neurons in the visual cortexOSB 

Which mechanisms are capable of generating a collective synchronous behav- 
ior? It is accepted that in order to observe a global coherent activity, oscillatory 
interacting units are required. The rhythmic activity of each element is provoked 
either by some internal processes or by external sources (external stimuli or forc- 
ing). The internal processes may have a physical or biochemical origin of great 
complexity, which probably is different for any of the systems considered above, 
but the essence of the phenomenon can be explained in terms of basic principles 
which allow to create a common framework. Within this environment, a simple and 
general mathematical description can be formulated. 

Suppose that the rhythmic activity of each element can be described in terms 
of a magnitude which evolves regularly in time. When such variable reaches a cer- 
tain threshold value the unit emits a pulse (action potential for neurons) which 
is transmitted to all the members of the population connected to it. Moreover, a 
resetting mechanism initialize the state of the unit that has fired. Then, a new 
cycle recommences. Essentially the behavior is analog to an oscillator. Assuming 
that the period of the rhythm is T, it is convenient to define the concept of phase, 
without loss of generality defined in [0, 1], which is a periodic measure of the elapsed 
time. In general, the phase is a nontrivial function of the state of the considered ele- 
ment. These features define the so-called relaxation or integrate-and-fire oscillators 
(IFO's). 

Essentially, the effect of the emitted pulse is to alter the current state of the 
neighbors by modifying their periods. This change is not uniform, sometimes the 
period is lengthened and sometimes it is shortened. The perturbation depends 
on the current state of the oscillator which receives the external impulse. The 
modification induced by such a pulse can also be studied in terms of a phase-shift, 
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what has been called phase response curve (PRC). Several different PRC's 
reported in the literature depending on the particular system under study, 
general, the analysis of the collective behavior of the system can be done in terms 
of a PRC if the phase-shift provoked by an impulse is independent of the number of 
impulses arriving within an intcrspike interval and if the arrival of such an impulse 
affects the period of the current interval but does not modify future intervals. This 
fact fostered a bunch of papers published in the seventies and the beginning of the 
eighties devoted to analyze the different patterns of entrainment between pairs of 
pacemakers. However, a systematic treatment of the whole population was still 
missing. 

There is another type of models where synchronization effects have been stud- 
ied extensively. In these models each member of the population is modeled as a 
nonlinear oscillator moving in a globally attracting limit cycle of constant ampli- 
tude. These are the phase, or limit cycle, oscillators. The units interact weakly to 
ensure that any perturbation does not move any of them away from the limit cycle. 
Then, only one degree of freedom is necessary to describe the dynamic evolution of 
the model. Perhaps, the best known example is the Kuramoto model@ where the 
phase (fi of each oscillator obeys the following Langevin equation 

N 

tpi =u>i + 22 Jij sin(ipj -ipi) + »?»(*) i (16) 
3=1 

where u>i is the intrinsic frequency of each member, and r]i denotes a gaussian white 
noise of zero mean and correlation given by 

< r)i(t)Vj(t') >= 2DS tJ 6(t - t'). (17) 

When there is no dispersion in the values of u>i this model can be transformed 
into the planar XY model, well known in statistical mechanics. However, in the 
context of oscillators u>i is picked up from a random distribution. Note that the 
main difference between these phase-coupled oscillators and the pulse-coupled IFO's 
explained before is that now the coupling in not pulsatile but it acts continuously 
in time, depending only on the phase difference between units. 

For ferromagnetic all-to-all interactions (Jij — J/N,J > 0), there is a critical 
value of the coupling strength J for which a phase transition occurs between a state 
where all the oscillators run incoherently with its own frequency to another state 
where a certain degree of synchronization appears spontaneouslv, The nature 



of such transition depends on the distribution of frequencies, EJtpon the features 
of the coupling,@0Q'lll and on the strength of the noise ZXeB In this review 
our interest will be focussed in the pulse-coupled model although under certain 
conditions it is possible to convert a model of IFQ's into a model of phase-coupled 
oscillators as has been shown by several authors E10O 

In general, three different levels of synchronization can be observed when con- 
sidering populations. The strongest implies that all the oscillators have the same 
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phase. When the members of a group do not fire in unison but keep a fixed phase 
difference between them, a "phase locking" regime is achieved. Finally, the weakest 
situation appears when several oscillators run with the same average frequency but 
without keeping any fixed phase difference. This is called "frequency locking". 

A particularly interesting task was developed by Winfree who was one of the first 
authors who considered_thc population as a whole and who gave mathematical tools 
to tackle this problem.L3 His ideas were picked up by other scientists who developed 
them more deeply or proposed alternative mathematical techniques to solve prob- 
lems which remained open for decades .l3 Mean- field models are sufficiently simple 
to be solved exactly and analyze under which conditions any of the three levels of 
synchronization defined above can be observed. 

3.1. Mean-field models 

We will consider systems formed by an assembly of N identical IFO's with all-to-all 
couplings. Each member of the population is described in terms of a phase variable 
ip and a state variable E which evolve in time according to the following equations 

t - 1 <** 

and 

HE 

-± = f(E i )+g{E i ,t), (19) 

where f(Ei) is the driving rate, which gives the natural evolution of the oscillator, 
and g(E{, t) is a function that accounts for the effect of the coupling between unit i 
and the rest. It may have a complicated dependence on time. In this description self- 
interaction is not considered but an explicit dependence of g(Ei,t) on the current 
state of the oscillator which receives the incoming input. The dynamics follows the 
essential trends mentioned in the introduction for relaxation oscillators. When an 
oscillator reaches the threshold it is reset [E — > 0, if — ► 0) and the rest of oscillators 
suffer a perturbation of their internal state. The features of this perturbation depend 
on g{Ei, t). This description is quite general and includes the majority of models of 
integrate-and-fire oscillators studied so far. Moreover, it allows a straightforward 
generalization to models with nonidentical periods. We can see in this dynamics a 
first contact point with models which display SOC behavior. 

A particularly interesting case was considered by Mirollp_and Strogatz (MS)Bj 
Their work was inspired in the analysis developed by PcskinEJ for the cardiac pace- 
maker, though the main outcomes are general enough to be applied to a wide range 
of models. MS consider systems whose intrinsic dynamics is given by a driving rate 
f(E) that satisfies 

f(E)>0 and /'(£) = fgl < , WE, (20) 

that is, the state variable E and the phase (p are related through a convex function 
E(ip), hereafter called the driving (see Fig. ||). 
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Fig. 5. Functional relation between the state of a unit E, and its phase ip, for different driving 
rates corresponding to a convex driving (f'(E) < 0), linear driving (f'(E) = 0), and concave 
driving (f'(E) > 0), respectively. 

As an example, for the Peskin modeliil 

f{E)= 1 {C-E), (21) 

where 7 is a constant that gives information about the convexity of the driving 
and C = 1 _ 1 e -- l ensures that the period is one when the threshold is one. This 
simplified model accounts for the variation of the membrane potential E (voltage 
difference between the internal and external part of a cell), as a consequence of the 
transport of Na, K and other ions across the membrane £hannels. The system is 
simply represented as an RC circuit, which implies 7 > 0e3e3 

On the other hand, in the MS model, the coupling between the oscillators is 
given by 

g{E i ,t)=eJ2Kt-t j ), (22) 

3 

where the sum involves all the oscillators which have fired at a given moment (tj 
denotes the time at which the j-th oscillator fires). An equivalent description can 
be made in terms of the state variable E 

g(E h t) (xe^2$( E c - Ej), (23) 
3 

where now the sum runs over all the oscillators and E c is the threshold. The 
coupling is instantaneous, excitatory (e > 0), and constant since there is no depen- 
dence on Ei , which means that when a cluster of m members reaches the threshold 
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simultaneously, the state of the other units is changed immediately according to 

E -> E + me. (24) 

One of the effects of such pulse-like interaction is that when one oscillator fires 
other members can also reach the threshold. These oscillators may in turn push 
other elements up triggering an avalanche. In this case all the units which have 
fired merge in a group which does not break up anymore. This effect is called 
absorption. The concept of absorption implies the existence of a refractory time, a 
period where the oscillator which has fired is not sensitive to new incoming inputs. 
It has a crucial effect when the interaction is modeled in terms of a delta function, 
as we will see later. Notice that a coupling given by (|2^) implies the existence of 
two time scales, a slow scale for the intrinsic dynamics of the units and a fast one, 
indeed infinitely rapid when compared with the other, for the interaction between 
IFO's. Once again we can establish a connection with models displaying SOC. In 
terms of this time scale separation we can say that the refractory time only acts in 
the fast time scale. 

MS were able to prove that for any e > the population of oscillators will be 
firing in synchrony in the stationary state for almost any initial condition (except 
in a set of Lebesgue measure zero) provided the driving is convex, i.e., when Eq. 
( p0| ) holds. In fact, this result can be derived even with more restrictive conditions 
given by 

g{Ei,t)=eS{t-tj). (25) 

According to this description the strength of the pulse emitted by a group of oscil- 
lators is independent of its size, that is, there is no additivity in the coupling and 
therefore m = 1 in Eq. (|24|). This result is very powerful because (|2^) does not 
favor synchronization, just the opposite. When additivity is taken into account, 
the process of entrainment is accelerated and perfect synchrony can be found under 
broader conditions than the given by the MS theorem.Eirc-3 

Note that there is no contradiction in these results because the theorem only 
gives sufficient but not necessary conditions to synchronization (Eq. (|20|)). Figure 
^ illustrates this situation. 

The MS model has been generalized to other situations. It has been^itudied 
in presence of weak noiseE3 and with a random distribution of thresholds. E3 Some 
studies have emphasized the necessity to overcome some unrealistic features of the 
model. For instance, the instantaneous character of the coupling does not take into 
account neither the finite time of the synaptic response nor the finite time associated 
with the transmission of information propagating along the axons (delays). Both 
are re 
richer. 



in models of spiking neurons and make the collective behavior much 
As an example, recent studies have shown that under cert ain con di- 
tions inhibitory couplings rather than excitatory lead to synchronization^ in 
agreement with some experimental observations.^ These results contrast with the 
MS model for inhibitory coupling (e < 0). A straightforward revision of their proof 



S8U 



Self- Organized Criticality and Synchronization in Lattice Models 21 



100 



80 



60 



to 



"i 1 1 1 1 1 1 1 r 



(a) 




20 40 60 80 100 120 140 160 180 200 



-t-a 




600 



Fig. 6. Time evolution of the size of the "avalanches" for the Peskin model with all-to-all coupling, 
e = 0.1, and TV = 100. (a) A convex driving, 7 = 0.1, leads to synchronization, and (b) with a 
concave driving, 7 = —0.1, the system settles in a periodic state which is very sensitive to initial 
conditions. In these simulations we have considered the extreme situation given byEq. (§|). Time 
is measured in number of avalanches. 
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shows that the sufficient conditions to find synchronization in this case are f(E) > 
but f'(E) > 0, VE. 

Another underlying hypothesis that restricts the range of interest of the MS 
model concerns the strength of the coupling e which does not depend on the state 
of the oscillator which receives the pulse. As we have mentioned above, in many 
studies of biological pacemakers it is assumed that the response of a cell due to 
an external stimulus is a phase shift (given by the PRC) which induces an energy- 
shift e(E) in Eqs. (|22]) or (pffi), hereafter called the energy response curve (ERC), 
that, in general, is a nonconstant function of E. It is possible to generalize the MS 
theorem by considering a general ERC as well as an arbitrary driving rate, with the 
only restriction of having f(E]> 0. Then, a sufficient condition to find a stable 

by§ 

< e'(E), WE. (26) 



synchronous regime is given bj 

f{E)e{E) 



f(E) 

This means that, in general, the knowledge of the intrinsic dynamics (i.e. the 
driving) of an individual oscillator is not enough to predict the final state of the 
population. Information about the response of a given unit in presence of a external 
stimulus is also required. 

3.2. Short range models 

The question whether the conclusions extracted from the MS theorem can also 
be applied to lattice models with short-range connectivity is of great importance. 
First, notice that in mean-field models synchronization emerges in an absorption 
process, where clusters of oscillators of increasing size merge with each other and 
never break up. However, in a lattice model, big assemblies of oscillators with the 
same phase, which eventually may break up, are generated by means of avalanches 
that start at a given point and propagate through the lattice. When they sweep the 
whole system we call them relaxation oscillations. Therefore, due to the different 
nature of both mechanisms one would expect that the conditions required to find 
synchronization/relaxation oscillations are different. In general, it is not easy to get 
analytical results for short-range models. Only for very special be andfor systems 
of low dimensionality (essentially ID) exact results have been derived.E3 Therefore, 
in the majority of situations the main outcomes rely on computer simulations. As 
we will see later, they show that under certain conditions mean-field arguments can 
be extrapolated to systems with finite connectivity. 



Let us consider a lattice model of IFO's with local rules given by (19) and ( |22| ) 
(but now the sum runs only over nearest neighbors) . An important point concerns 
the way in which an avalanche is propagated through the lattice. Suppose that a 
given unit, the seed, has fired and as a consequence some of its neighbors reach the 
threshold. In turn, these oscillators will fire, interacting with their neighbors being 
the seed among them. The question is what happens with the seed. According to 
the strategy applied in the mean-field version of the MS model such unit should 
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remain at the reset point during the whole avalanche. This is important because 
such effect implies the existence of a refractory time which persists until all the 
oscillators are below the threshold. Furthermore, all the members which have fired 
will have zero phase when the avalanche is over. We can describe this refractory 
time by an ERC that verifies e(0) = 0. Under such conditions computer simulations 
show that relaxation oscillations appear when f'(E) < and_£ is a definite positive 
constant, in agreement with a conjecture proposed by MSEJ For a more general 
function e(E) (but e(0) = 0), it has been showno that (|26| ) is a sufficient condition 
to observe relaxation oscillations. Moreover, in this case relaxation oscillations lead 
to a perfect synchronous regime. c Figure |?] gives numerical evidence of this fact. 
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Fig. 7. Number of avalanches (filled) and time in period units (hollow) required for a 32 X 32 
lattice of Peskin's oscillators with open be to observe complete synchronization as a function of 
7 in a log — log scale. The straight line shows the 7 behavior. The value of e is constant and 
equal to 0.1. Similar results are obtained for periodic be and other driving rates. 



The most appealing consequence of the previous results is that for coupled dy- 
namical systems like those considered before, for open or periodic be, mean-field 
results seem to be still valid and represent sufficient conditions to ensure the exis- 
tence of the synchronized regime, provided one considers the refractory time. 

When the driving is linear, the phase and the energy become identical. In this 
case, for a constant ERC, it is very simple to show that starting from a random 
initial configuration the system evolves towards a periodic state formed by clusters 
of oscillators with the same phase. The phase difference, after a cycle, remains 
constant in time once all the units have fired, and consequently the system does 

c Note that for open be and e(0) 7^ 0, after a relaxation oscillation of the size of the system not 
all the units have the same phase. For instance, the cells located at the boundaries remain at a 
different state with respect to those located in the bulk because the number of neighbors in each 
case is different. 
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not synchronize in agreement with Fig. Q]. On the other hand, a concave driving 
may give rise to the appearance of stable complex phase- locked spatio-temporal 
structuresEJEj'Ell as, for instance, a chessboard- like configuration. This sort of "an- 
tisynchronized" states are exclusive of the short-range models since there is no 
counterpart in the mean-field case. 

On the other hand, we could suppose in the previous discussion that the seed 
was updated as any other neighboring cell, that is, e(0) ^ 0, so the refractory time is 
removed. In that case the coupling between neighbors is governed, when > E c , 

by 

These rules are not new. They were introduced in a stick-slip model by Feder and 
Feder (FF)l3 in the context of earthquake dynamics, with the restriction of having 
a constant driving rate and open be. 

As in all the models of IFO's discussed so far, since the relaxing site is reset 
to zero and a fixed amount is transferred to the neighbors, the FF model does not 
have any memory of the initial configuration. After an avalanche, all the sites that 
have fired have an energy which is a multiple of e. When the threshold E c and e 
are commensurable, that is E c /e is an integer M, the system eventually reaches a 
state with M distinct values of the dynamical variables. Hence, after an avalanche, 
Ei_j = ne, n = . . . M — 1, and since the driving is uniform, Ei_j — ne, n = 1 . . . M 
when initiating an avalanche. When E c is incommensurable with e, that is, when 
E c /e is not an integer, the system is not able to reach a "partially synchronized" 
state. Figure |8|(a) displays the average number of totally synchronized sites as a 
function of time for E c /e — 4, 100/22, and 5. When E c is commensurable with e, 
the time the system needs to settle down into such a state increases as e decreases. 
Note that the case of M — 4, was originally studied by Feder and Feder .ta Their 
results are quite different. This will be explained in detail in the next section. Figure 
||(b) displays the corresponding distributions of avalanche sizes. When M = 4 and 
5 only large avalanches are seen. For M = 100/22, avalanches of all sizes occur, but 
they are not power-law distributed. 

4. Synchronization and Self-Organized Criticality 

Up to now, we have seen the strong similarities between some models displaying 
SOC and lattice models of IFO's showing different levels of synchronization. In 
both types of systems, each unit has an intrinsic dynamics leading to a threshold 
value. When this threshold is reached, the unit interacts with its neighbors and it 
is reset. This process is developed in a time scale much faster than that associated 
to the natural driving of each unit. The features of the couplings, the natural 
dynamics, and the reset mechanism define each particular situation. In view of 
this analogy it is reasonable to think that a general framework may be developed. 
We have shown that by an appropriate choice of the relevant parameters of these 
models a wide spectrum of collective behaviors can be observed. However, up to 
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Fig. 8. (a) The average number of synchronized sites in the FF model with e = 0.20, 0.22, and 
0.25, respectively. E c = 1 and system size L = 50. (b) The corresponding avalanche distribution 
when the system has reached the final partially-synchronized state. Time is measured in number 
of avalanches. 
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now our attention has been focusscd on systems which display either SOC or any 
sort of phase locking. The goal of this final section is to give a step forward by 
showing that there is an interesting group of models, with quite different origins, 
which display both behaviors. 



4.1. Integrate- and- fire oscillators 

Let us start the discussion with a model studied in the previous section. We have 
seen that the FF model tends to develop, in the stationary state, either assemblies 
of cells with the same phase (when the ratio E c /e is an integer) in which case the 
big avalanches take place because a lot of units reach the threshold simultaneously 
or a more complex behavior (when E c /e is not an integer) where avalanches of 
all sizes take place, hat they are not power-law distributed. These behaviors are 
not robust to noiseEll In principle, the effect of the noise is to prevent that two 
sites become critical simultaneously due to the uniform driving. The nature of such 
noise can be based on the existence of fluctuations: on the dynamic variables, on 
the thresholds, or on the reset energy.Eilo This source of dynamical noise has to be 
distinguished from a quenched source of noise, as for instance, a random distribution 
of frequencies. 



UK) 



The dramatic effect of noise on the features of the stationary state of the model 
is observed in Fig. ||. When altering the relaxation rule (|27]) slig rUy by adding a 
small random number in the range [0, 0.001] to the reset oscillator the system no 
longer settles into a state with a few groups of oscillators with the same phase, but 
it goes towards a new state in which the avalanche distribution follows a power-law 
decay characteristic of SOC. 

Keeping in mind the results given in the previous section for the MS model, it is 
not difficult to imagine that by increasing the convexity of the driving (in the original 
FF model it is linear) SOC fades out, even in the presence of a dynamical noise, d 
and for a certain degree of convexity a spatial structure formed by large groups 
of oscillators with the same phase appears spontaneously.E2i In fact, the features 
of the model allow to derive analytically an explicit condition, written in terms of 
a critical convexity j c , which ensures that a state characterized by the existence 
of relaxation oscillations of the size of the system is a fixed point. However, this 
argument cannot ensure that such fixed point is an attractor. We will show, by 
numerical simulations, that it has a large basin of attraction. 

Let us assume that there is an avalanche sweeping the whole lattice (for an arbi- 
trary geometry). The idea is to derive a condition that could ensure that the next 
event will reach the whole population as well. Let n denote the lattice coordination 
number and m the number of neighbors of the most unfavorable site (the cell or 
group of cells with the smallest number of neighbors). These sites have zero energy 



d Corral et aiJiiiU introduced the noise by choosing at random one of the cells that become critical 
simultaneously due to the driving and making it the starting point of the avalanche, while the 
energy of the rest of critical cells was decreased by a small random amount. However, different 
sources of dynamical noise give similar results. 
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when the avalanche is over, and in order to fire in the next avalanche, they have to 



satisfy££il (for excitatory coupling), 

E(l - <f(ne)) +me > E c . (28) 
In the case of the Peskin model the relations between the phase and the energy are 

E(<p) = C (1 - e" w ) (29) 

and 

= i ha (30) 
where C was already defined in (|l|). A simple calculation yields, when E c = 1, 

e < C(m '" ) + n (31) 
mn 

which is the relation between e and 7 that must be satisfied for model ( pl| ) with in- 
teraction rules (|27]) in order to show relaxation oscillations involving all the cells. In 
the particular case of a regular 2D square lattice with open be, the most unfavorable 
sites are the corners. Hence 

e < 2 ~^- (32) 

This relation has been checked through simulations finding an excellent agreement as 
can be seen in Fig. |l^. In region A, the condition ([32]) is satisfied and therefore, the 
system settles in a state where a macroscopic group of oscillators whose size scales 
with the size of the lattice remain at the same phase. In region B, the condition is 
not satisfied and relaxation oscillations involving all the cells are not stable. The 
systems goes towards an attractor of difficult classification whose characteristics 
will be discussed later. In the limit of 7 = we recover the features of the noisy FF 
model. Notice that a convex driving is not enough to find relaxation oscillations. 
It is necessary to go beyond 7 C . In fact, 7 C tends to In 2 in the limit of e — » 0. 
When considering periodic be in a 2D square lattice, Eq. ( |3l| ) reads 

e < 0.25 (33) 

thus favoring the appearance of relaxation oscillations in the sense discussed above: 
it is a stable fixed point but, in principle, we cannot ensure it is an attractor for the 
dynamics of the system. 6 Indeed, numerical simulations show that this synchronized 
state can be an attractor of the dynamics in some cases. This model has been also 
analyzed by Hopfield and Herzt£3'll2^l in the context of computation by networks of 
integiale-and-fire neurons. These authors study a general model with interaction 
ruleai^ 

Erin > E nn -\- (34) 



e For e > 0.25 the model with periodic be has no dissipation and finally enters into a never-ending 
avalanche. 
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Fig. 10. Schematic phase diagram in terms of 7 (the convexity of the driving) and e (the strength 
of the coupling) for model (gfl) and units driven by (pi]). The symbols correspond to the phase 
transition observed in the simulations and the error bars are given by the standard deviation over 
ten measures. The solid line corresponds to the analytical result (|32[) . 




and a linear driving. This model is equivalent to the FF model in the T = limit 
and to the uniformly driven BTW when T = 1. As soon as every element has 
fired once, the system reaches a limit cycle in which each oscillator fires exactly 
once during one period of the oscillation, no matter the value of T. For the BTW 
model it is possible to construct a Lyapunov function that show this convergence. 
However, for T < 1, the volume of the attractor is greatly reduced. From numerical 
simulations on 2D square lattices they construct a diagram in terms of the initial 
distribution of energies and the coupling strength e. The behavior ranges from a 
trivial fully synchronized state, when the initial distribution of the oscillators is 
very narrow, to an exponential distribution of events, when the initial distribution 
is broad enough. In the border line between both regions a power-law distribution of 
events is reported. The same authoraI23 have also studied the limiting cases of the 
BTW and FF models with units being driven nonlinearly. Within the computational 
aspects of this work they perform simulations from which one can conclude that all 
the models exhibit rapid convergence to periodic solutions with locally synchronized 
clusters. However, in the large time scale the behavior is quite different. In this 
case, the BTW interaction rules lead to a global synchronization whereas FF's show 
a partially synchronized behavior. The existence of this nonglobal synchronized 
periodic state for the FF model confirms that Eq. ( |33| ) gives the existence of a 
synchronized state (in the sense of a simultaneous firing of all the units) that, 
however, it is not necessarily a global attractor for the dynamics. These authors 
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give conditions which show that if the strength of the coupling is strong enough, a 
simple spatially nonuniform periodic solution is also possible. 

4.2. Spring-block models 

Now, we are going to analyze the effect of a nonlinear convex driving in deter- 
ministic spring- block models, such as those discuss ed- in -the context of earthquake 
dynamics, which display SOC with a linear driving.ll3e3't£il As in models of IFO's, 
these nonlinearities tend to break SOC in favor of more complex dynamical states. 
However, there are several differences that must be remarked. In the OFC model (7) 
two cells very seldom reach the threshold simultaneously and trigger an avalanche 
in two different points of the lattice. There is no need to introduce a dynamical 
noise in contrast with the FF model where the existence of the noise is required 
to observe SOC. Therefore, concerning the SOC behavior, one might conclude that 
the existence of some kind of noise, either in the way the system is driven or in 



.bora 



105 



iave reported 
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i£3 However, 



the initial conditions, is necessary. On the other hand, some aut 
power-law distributions for completely deterministic automata, 
for all these deterministic models there is a single site with a special behavior; this 
inhomogeneity, either in the initial conditions or in the external loading, is what 
gives rise in these cases to the complex behavior characteristic of SOC. 

For the OFC model it is also possible to deduce analytical results giving in- 
formation about the conditions required to observe a transition between a regime 
where relaxation oscillations rule the long time behavior and other dynamical states. 
However, the situation here is more difficult and two equations are necessary. One 
ensures that the avalanches will reach all the cells of the lattice, and the other is 
based on the assumption that the configuration after a big event always will be the 
same. The main difference with respect to the FF model is that we have to replace 
e by an effective value s since the energy of a given site can be larger than E c = 1 
when the avalanche propagates through the lattice/ Within a mean-field approx- 
imation we will assume that this energy is the same for all cells, except the seed 
(the unit that triggers the avalanche), its neighbors and the boundaries. Assuming 
that any site fires at most once during an avalanche 3 the seed has a phase <p(4e) 
when the avalanche has finished, and to fire again it has to increment by an amount 
1 — </?(4e). Since its neighbors are at tp(3e) the condition to repeat permanently this 
situation is 

4e = 4e [E (y>(3e) + 1 - <p(4e)) + e] . (35) 

This equality, in addition with (^Tj) , with e replaced by e, gives the following rela- 
tionship between s and 7 



e < — 
~ 16 



VC 2 - 52C + 164 - (C + 6)1 . (36) 



^Actually, for a regular 2D square lattice e is bounded between e and 2e since the energy at any 
site cannot exceed the value of 2. 1—1 
9 This can only be ensured when e < 0,22.r 9 l 
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When this condition is fulfilled, only relaxation oscillations of the size of the system 
can survive in the stationary state. The curve corresponding to the equality has 
been plotted in Fig. [ll] (solid line). It has been corroborated through simulations 
performed in a system of size 64 x 64 with a fixed value of 7 and increasing e from 
below. For 10 different realizations of the initial random conditions we look at the 
final state after 3 • 10 6 avalanches. The average value of e where the transition is 
observed corresponds to the circles in Fig. |ll| and the error bars correspond to the 
standard deviation. 
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Fig. 11. Schematic phase diagram in terms of 7 (the convexity of the driving) and e (the strength of 
the coupling) for the OFC interaction rules with a driving given by (pl[). The symbols correspond 
to the phase transitions observed in the simulations. For the B-C transition the error bars denote 
that above them we have always found SOC while below them there is no power-law behavior. 
For the A-B transition the error bars are given by the standard deviation over ten measures. The 
solid line corresponds to the analytical result (tSq) and the dashed line is an exponential fit to the 
numerical data. 



Furthermore, to complement and illustrate the features of the dynamical states 
observed in the phase diagram of the model, we have plotted in Fig. 



12 the time 



evolution of the avalanche size as well as its return map, for a particular run. Thus, 
we can clearly see that, in region A (Fig. |l^(a)) large avalanches appear frequently 
until an avalanche sweeps the whole system and this state is maintained. On the 
other hand, the return map shows the evolution of the system at the fixed point. Let 
us remind that in the theoretical approach we made the assumption that a whole 
system avalanche does exist. In the simulations one observes that this hypothesis is 
indeed true. The synchronization of the cells is not global but the number of cells 
with the same phase is of the order of the size of the system. 

Slightly above the curve given by Eq. (|36| ) an avalanche sweeping all the lattice 
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cannot repeat any more since the next one will be unable to reach some of the cells 
in the boundaries and these cells will be the starting point of future avalanches 
(region B). This fact gives rise to a periodic behavior with a discrete distribution 
of a few avalanche sizes that can eventually be broken only by the effect of the 
dynamical noise. This is indeed what we have observed in the simulations; only 
some values of avalanche sizes are present. The exact location is very sensitive to 
the initial random configuration but they are always of order 1, L, and L 2 . In 
Fig. |l2|(b) we show the time evolution of the avalanche size in a typical run in this 
region after a transient time. Here we can notice the appearance of the periodic 
state and a small number of points in the return map. 

By increasing e we observe that the distribution of avalanche sizes, P(s), be- 
comes continuous in the sense that all bins of exponentially increasing width have 
at least one event, up to some system size dependent cutoff. The transition from 
a sparse distribution to a "continuous" one would require a careful investigation. 
When P(s) has this continuous appearance there are peaks at some characteristic 
lengths, roughly at i, 2L, 3i, . . ., as can be seen in Fig. |l3], and no finite-size 
scaling is possible. 

By increasing e again we observe that the intensity of the peaks decreases. Up to 
the system sizes and number of avalanches we have studied, there exists a transition 
to a stationary state with a power law, followed by an exponential decay of P(s). 
We have identified the transition from regions B to C (squares in Fig. |ll|) by fixing 
7 and increasing e up to the appearance of a finite-size scaling in the distribution 
of avalanche sizes. The dashed line is an exponential fit of the numerical results 
that we extrapolate to the linear driving case (7 = 0) and for large values of the 
convexity of the driving. 

The criticality of the SOC state stems from the lack of characteristic time and 
length scales and for the apparent unpredictability of the size of the next event. This 
can be visualized in Fig. [l2|(c) where no simple correlation can be extracted from 
the time sequence. Just at the conservation line, s = 0.25, where the interaction 
between units is stronger, we find SOC for a wide range of values of 7 with a 
power-law decay over more than three orders of magnitude and the corresponding 
finite-size scaling as shown in Fig. 

Up to now we have studied a population of identical oscillators, with intrinsic 
period T = 1, but in a real system it is usual to find it randomly distributed. 
We have consider a uniform distribution, centered at T = lJl23 For model (|9|) 
and a convex driving the results remain qualitatively the same for a width in the 
distribution as large as 0.1T. However, if the width is equal to the mean period 
T the region where synchronization is observed almost disappears as well as SOC 
which is constrained to a narrow region near the line e = 0.25. The effect of a 
quenched distribution of thresholds has been also analyzed by Janosi and Kertesz 



(uniform distribution)E2a and by Torvund and Fr0yland (gaussian distribution) 



for the linearly driven OFC model, where it has been shown that the distribution 
of avalanche sizes has an exponential decay for a sufficiently large spread. 
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Fig. 12. Left column: Time evolution of the avalanche size. Right column: return map of the 
avalanche size. Top) in region A of Fig. [ll] for L = 8, middle) for region B (close to region A), 
and bottom) for region C. Each time step corresponds to an avalanche. 
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Fig. 13. Distribution of avalanche sizes in region B (close to region C). 
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A slightly different model with a uniform driving was analyzed by Socolar, Grin- 



stein, and Jayaprakash.LL2-3 In this case when a site reaches a critical value the in- 
teraction rules are 

Enn * A-£/ nn -(- £Ei j /o>7\ 

25„ - 0. (37) 

In this model, A ^ 1 is a simple way to remove degeneracy, though it does not 
model the elastic nonlinearities of the blocks and springs model in detail. Rules 
(|37|) can be seen as a kind of a state dependent transfer (ERC), which can be 



transformed into a uniform 
appropriate transformation 



')!!.')] 



jfer model with a nonlinear driving by means of the 
When open be are considered, A = 1 corresponds 
to the original OFC model and in an interval around this value there are indications 
of SOC behavior. Nevertheless, for A < 1, the authors suggest the observed behavior 
might correspond to a periodic state with either a long transient or a long period. 
For A > 1, the power law decay in the avalanche size distribution observed is related 
to apparently chaotic oscillations in the average energy 

<£>=-^X>,i. (38) 

It is suggested that the boundary sites, which have a different cycle than interior 
sites, might be the source of the criticality observed. This issue has been analyzed 
in detail by Middleton and TangEH in a simple directed model; they show that the 
slower growth rate of the boundary sites favors the creation of synchronized clusters 
of all sizes. At a fixed time, the number of clusters of size d has a power-law tail 
that is independent of e, n(d) ~ d~ a . An avalanche can then eventually cross a 
cluster boundary whenever both clusters have close enough values of the energy, 
and this makes the avalanche size distribution to be also a power law. 

As we noticed for the FF model, also for the OFC model a very important role 
is played by the be. In particular, periodic be break completely the SOC behav- 
ior and can p rovide s ome hints on its origin in the nonconservative continuo usly 



T i°n ii 1' 



104 



driven models .c-3liii3l£j For instance, Socolar, Grinstein, and Jayaprakash obtain 
different behaviors, depending on the values of A in (|37|). Thus, for A < 1, the 
system always settles in a periodic state, with transients that become longer as A 
approaches one. On the other hand, for A > 1 the most common states are peri- 
odic with a cycle consisting of single avalanches that sweep the whole lattice. A 
simple picture can be obtained by analyzing the return map of a two-site model: 
for A > 1 both sites fire in the same avalanche, whereas for A < 1 there is a stable 
fixed point. These results are in agreement with the onesobtained for a nonlinear 
driving by means of the above mentioned transformation.^ Finally, the special case 
of A = 1, shows a completely different behavior from what is expected from the 
limits of the previous cases: the system rapidly settles into a periodic state with 
avalanches of size one which are marginally stable as those described previously for 



the FF model Jl23'E23 Grassberger,EJ by means of very large-scale simulations, real- 



ized that there exists a close relation between temporal and spatial ordering; while 
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time ordering is associated with the periodic behavior and becomes stronger as £ 
decreases, spatial order is related to the synchronization of neighboring sites which, 
in principle, seems to be favored by a large interaction between units. Concerning 
the SOC behavior, one can conclude that open boundaries create an inhomogeneity 
which prevents time orderin g h ut spatial local ordering is still maintained. Along 
this line Middleton and Tango have observed that by including some source of dy- 
namical noise the system with periodic be is periodic in time, after a brief transient 
time, and that these states are neutrally stable. 

A different way to introduce the inhomogeneity th at is responsible for the SOC 
behavior has been followed by Torvund and Fr0y land They have considered a 
single site to have a threshold larger than the threshold of the remaining sites, 
in a system with periodic boundary conditions. These authors have found some 
evidence for a SOC behavior when the larger threshold is above a critical value. 
Nevertheless, also a non-simple periodic behavior can be observed. This means 
that a more careful analysis would have to be performed along this line. 



4.3. Other models 

Forest- fire_. models are other systems where typical SOC behavior has been 
observed.^ The model is defined in a £>-dimensional lattice (usually D = 2). Each 
site of the lattice is either occupied by a tree or empty. Trees grow according to 
a certain function that depends on the age of the empty site. A tree burns if it 
is struck by lightning with a small probability /. In this case, the fire propagates 
through the neighbors burning the whole cluster that immediately becomes empty. 

In this system the features of the stationary state depend strongly on the specific 
form of the tree growth. Drosse£L3 has shown analytically (for a ID system) and 
through simulations that if the life-time distribution of empty sites has a long-time 
tail a SOC state with non-universal exponents is observed. On the other hand, 
for deterministic tree growth a coherent temporal activity between trees in the 
same cluster comes up. Therefore, we see again the relevance of the driving in the 
collective properties of the model. 

Inspired by a forest-fire model, Clar, Drossel, and Schwablli^l have studied a 
nonequilibrium percolation system which displays several behaviors ranging from 
SOC to clusterization. The authors consider a population of objects (trees, animals 
or others) distributed randomly over a lattice. The density of occupied sites is p. 
Then, a given site is chosen. If it is not empty an "explosion" occurs affecting 
all the neighbors within the same cluster, whose spatial position on the lattice is 
redistributed (a group of animals dispersed by the action of a predator could be 
another example) . In these terms the density of the system is a conserved quantity. 

In the stationary state the model shows different regimes depending on the 
density of occupied sites of the lattice. For diluted networks (p < p\ = 0.41) only 
small clusters of units are developed. As a consequence, the average number of 
units affected by an explosion and the size of the largest cluster s max is small. 
For pi < p < p2 = 0.435 the size of the explosions diverge but more slowly than 
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the size of the system. They show that the relevant quantities diverge as a power 
law, with an exponent that depends on the density p. This critical behavior is 
characteristic of SOC. Finally, for larger values of the density there is a region where 
the explosions have infinite size (they scale with the size of the system) which is 
repeated periodically in time. 

Related to this model we can notice the work of Csilling et alr!3 about dynamics 
of populations. Each site on a lattice denotes a population that evolves according 
to 

N id {t + 1) = XN itj (t) [1 + aN tJ (t)}- p (39) 

until a critical population density is reached (overcrowding) . In this case a dispersal 
movement (migration) is triggered through nearest neighbor local habitats 

N nn ► N nn + A ^ (4f)l 
N ii:j ^N sc . K ' 

Notice that these rules are identical to OFC rules with N sc — and A = 4e. This 
process may lead to a migration avalanche, typical of the lattice models discussed 
so far. Without interaction, a simple site behaves chaotically. However, by increas- 
ing the interaction it reaches a noisy fixed point. Concerning the lattice average, 
by lowering the threshold level, the collective behavior becomes more pronounced, 
reflected by the appearance of discrete frequencies on the power spectra; more- 
over, the time evolution either becomes strictly periodic or it reaches a stable fixed 
point. Finally, the distribution of avalanche sizes is computed; whereas for a weak 
interaction the distribution is exponential, for a stronger interaction a power law, 
characteristic of SOC, was observed. 



5. Conclusions 

In this paper we have reviewed the collective behavior of a large group of low di- 
mensional systems characterized by two essential features: a dynamical process to 
update the state of each member of the system and a threshold condition which 
defines an interactive process between members of the population. The coupling 
between units is defined in terms of very simple local rules. Within this broad frame- 
work we have considered different models ranging from cellular automata where the 
state of each unit, described in terms of discrete variables, is updated through a 
random process, to coupled dynamical systems where the internal state is governed 
by differential equations. In spite of the simplicity of the considered models, the 
majority displays a complex behavior. In some cases they settle down in stationary 
states without characteristic time or length scales typical of self-organized critical- 
ity (SOC). In other cases, the attractors are synchronized states with a partial or 
complete coherence in the temporal activity of the members of a given population. 
Most interestingly, a few of them display both attractors (and others) depending 
on the value of the parameters which define the model. 

There are several reasons which foster the study of these models. Perhaps, the 
most relevant is to understand the origin of SOC and identify the underlying mech- 
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anisms. Unfortunately, there is no clear answer to this question and the problem 
remains open. SOC behavior is very sensitive to the action of external factors, it 
can be destroyed (and created) in many different ways and, a priori, it is difficult to 
find the proper combination of ingredients that may make a system to display it. It 
is well known that boundary conditions play a relevant role. Big events (avalanches) 
are usually generated near the boundaries and propagated towards the bulk, but 
only under suitable conditions which, in general, imply open boundaries. But apart 
from this important factor, there are other components which also have a crucial 
effect on the features of the stationary state. For some cellular automata such as 
the BTW model conservation is the key. For the Zhang model cither conserva- 
tion or a very tiny perturbation of the internal state of the cells. Some continuous 
driven models display SOC even without conservation such as the OFC model but 
additionally, in others such as the FF model a dynamical noise is required. In the 
last two models the role of memory effects (a trace of the initial state) can give the 
answer. 

We have seen that in the continuous driven models the features of the driv- 
ing are also relevant. The original OFC and the noisy FF models display SOC 
with constant or quasiconstant driving rate. However, a non-constant driving rate, 
which implies a convex (concave) relation between the state of each unit and time, 
leads to new behaviors. In some cases the duality convexity /concavity lead to full 
synchronization/anti-synchronization (phase locked state where neighbors remain 
at a fixed phase difference). Then SOC might be understood as an intermediate 
behavior, balanced between these two extreme situations. In other situations, local 
synchronization can be responsible for a global SOC behavior. For these reasons, 
the analysis of systems displaying a rich variety of attractor such as those mentioned 
in the last section of the paper is of great interest. 

Another important feature of the systems considered in this review is the sep- 
aration of time scales associated to the natural dynamics of each unit and the 
interaction between them. This separation, quite natural in some situations, seems 
to be mandatory to observe SOC behavior. If this is true, systems where the trans- 
mission of information between neighbors could take a finite time might not display 
SOC. It would be interesting to pay more attention at this point in the future and 
see whether it is possible to observe power-law behavior in systems where both time 
scales are overlapped. 

Finally, there are other aspects of interest not sufficiently analyzed so far. The 
effect of different sources of noise (for instance, the coupling with a thermal bath) , 
the role played by inhibitory rather than excitatory couplings, a careful analysis of 
the role of memory effects are some of the topics that deserve further research in 
the future. 
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